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Tunnel transport processes are considered in a square lattice of metallic nanogranules embedded 
into insulating host to model tunnel conduction in real metal/insulator granular layers. Based on a 
simple model with three possible charging states (±, or 0) of a granule and three kinetic processes 
(creation or recombination of a ± pair, and charge transfer) between neighbor granules, the mean- 
field kinetic theory is developed. It describes the interplay between charging energy and temperature 
and between the applied electric field and the Coulomb fields by the non-compensated charge density. 
The resulting charge and current distributions are found to be essentially different in the free area 
(FA), between the metallic contacts, or in the contact areas (CA), beneath those contacts. Thus, 
the steady state dc transport is only compatible with zero charge density and ohmic resistivity in 
FA, but charge accumulation and non-ohmic behavior are necessary for conduction over CA. The 
approximate analytic solutions are obtained for characteristic regimes (low or high charge density) 
of such conduction. The comparison is done with the measurement data on tunnel transport in 
related experimental systems. 

PACS numbers: 73.40.Gk, 73.50.-h, 73.61.-r 



I. INTRODUCTION 

More than three decades have passed since the pioneer- 
ing studies by Abeles and co-workers [1, 2] that triggered 
a huge research effort in granular thin films. Actually, 
nanostructured granular films are of a considerable in- 
terest for modern technology due to their peculiar physi- 
cal properties, like giant magnetoresistance [3], Coulomb 
blockade [4, 5], or high density magnetic memory [6], im- 
possible for continuous materials. 

However, a number of related physical mechanisms 
still needs better understanding, in particular, trans- 
port phenomena in these films are still a great challenge 
and presently various works are addressing such prob- 
lem [7-9]. The main reason is that granular systems re- 
veals certain characteristics which cannot be obtained 
neither in the classical conduction regime (in metallic, 
electrolyte, or gas discharge conduction) nor in the hop- 
ping regime (in doped semiconductors or in common tun- 
nel junctions). Their specifics is mainly determined by 
the drastic difference between the characteristic time of 
an individual tunneling event (~ h/ep ~ 10~ 15 s) and 
the interval between such events on the same granule 
~ e/(jd 2 ) ~ 10~ 3 s, at typical current density j ~ 10~ 3 
A/cm 2 and granule diameter d ~ 5.0 nm. Other impor- 
tant moments are the sizeable Coulomb charging energy 
E c ~ e 2 /(e c gd) (typically ~ 10 meV) and the fact that 
the tunneling rates across the layer may be even sev- 
eral orders of magnitude slower than along it. The in- 
terplay of all these factors leads to unusual macroscopic 
effects, including a peculiar slow relaxation of electric 



charge discovered in experiments on tunnel conduction 
through granular layers and granular films [11, 12]. 

For theoretical description of transport processes in 
granular layers (and multilayers) we develop an exten- 
sion of the classical Sheng- Abeles model for a single layer 
of identical spherical particles located in sites of a sim- 
ple square lattice, with three possible charging states (±, 
or 0) of a granule and three kinetic processes: creation 
of a ± pair (the only process included in the original 
Sheng-Abeles treatment) on neighbor granules, recom- 
bination of such a pair, and charge translation from a 
charged to neighbor neutral granule. Even this rather 
simple model, neglecting the effects of disorder within a 
layer and of multilayered structure, reveals a variety of 
possible kinetic and thermodynamical regimes, well re- 
sembling those observed experimentally. 

The detailed formulation of the model, its basic param- 
eters, and its mean- field continuum version are given in 
Sec. II. Next in Sec. Ill we calculate the mean values of 
occupation numbers of each charging state under steady 
state conditions, including the simplest equilibrium situ- 
ation (no applied fields), in function of temperature. The 
analysis of current density and related kinetic equation 
in the out-of-equilibrium case is developed in Sec. IV, 
where also its simple, ohmic solution is discussed for the 
FA part of the system. The most non-trivial regimes 
are found for the CA part, as described in Sec. V for 
steady state conduction with charge accumulation and 
non-ohmic behavior. The general integration scheme for 
non-linear differential equation, corresponding to steady 
states in FA and CA, and particular approximations lead- 
ing to their analytic solutions are dropped into Appendix. 
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FIG. 1: Square lattice of metallic granules in the insulating 
matrix. 



II. CHARGING STATES AND KINETIC 
PROCESSES 

We consider a system of identical spherical metallic 
nanogranules of diameter d, located in sites of simple 
square lattice of period a within a layer of thickness b ~ a 
of insulating host with a dielectric constant e (Fig. 1). In 
the charge transfer processes, each granule can bear dif- 
ferent numbers a of electrons in excess (or deficit) to the 
constant number of positive ions and the resulting excess 
charge ere defines a Coulomb charging energy ~ <r 2 E c . 
At not too high temperatures, k-sT < E c , the consid- 
eration can be limited only to the ground neutral state 
cr = and single charged states a = ±1. Actually, for low 
metal contents (well separated, small grains) , E c reaches 
~ 10 — 30 meV, so this approach can be reasonable even 
above room temperature. For a three-dimensional (3D) 
granular array, E c was defined in the classic paper by 
Sheng and Abeles [1], under the assumption of a constant 
ratio between the mean spacing s and granule diameter d, 
in the form E c = e 2 / (s / d) / (ed) , where the dimensionless 
function /(z) = l/(l+l/2z). Otherwise, the complete di- 
electric response of 3D insulating host with the dielectric 
constant e and metallic particles with the volume frac- 
tion / < 1 and diverging dielectric constant e m —¥ oo can 
be characterized by the effective value e e // = e/(l — /). 

For the planar lattice of granules, the analogous effec- 
tive constant can be estimated, summing the own energy 
e 2 /(ed) of a charged granule at the n = site and the 
energy of its interaction with electric dipolar moments 
« (e/e e ff)(d/2n) 3 n, induced by the Coulomb field from 
this charge (in macroscopic dielectric approximation) on 
all the granules at the sites n = a(ni,ri2): 



E « = ~d 



1 



'eff 



e ef fd 



(1) 



Here the constant a = jYl 
suiting e e ff -- 



n=£0 ' 



« 5.78, and the re- 
■ + y/e 2 + ea(d/a) i /2 > e. However, 

Eq. 1 may considerably underestimate the most impor- 
tant screening from nearest neighbor granules at d ~ a, 
and in what follows we generally characterize the com- 
posite of insulating matrix and metallic granules by a 



certain 



■eff 



e 2 /dE c > e. 



Following the approach proposed earlier [11], we clas- 
sify the microscopic states of our system, attributing the 



charging variable er n with values ±1 or to each site 
n and then considering three types of kinetic processes 
between two neighbor granules n and n + A (Fig. 2): 

1. Electron hopping from neutral n to neutral n + A, 
creating a pair of oppositely charged granules: 
(cr n = 0, cr n+ A = 0) (er„ = +1, <7 n +A = -1), 
only this process was included in the Sheng and 
Abeles' theory; 

2. Hopping of an extra electron or hole from n to neu- 
tral n + A, that is the charge transfer: (er n = 

±1) C„+A = 0) -» (er„ = 0, cr n+A = ±1); 

3. Recombination of a electron-hole pair, the inverse 
to the process 1.: (a n = +l,o p n +A = ~ 1) ~~ ^ (°n — 
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FIG. 2: Kinetic processes in a granular layer. 



Note that all the processes 1) to 3) are conserving the 
total system charge Q = ^ n er n , hence the possibility for 
charge accumulation or relaxation only appears due to 
the current leads. A typical configuration for current-in- 
plane (CIP) tunneling conduction includes two macro- 
scopic metallic electrodes on top of the granular layer, 
forming contact areas (CA) where the current is being 
distributed from the electrodes into granules, through an 
insulating spacer of thickness 6', and a free area (FA) 
where the current propagates over the distance I between 
the contacts (Fig. 3). To begin with, we consider a sim- 
pler case of FA while the specific analysis for CA with an 
account for screening effects by metallic contacts will be 
given later in Sec. V. 




CA FA CA 

FIG. 3: CIP conduction geometry. 
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The respective transition rates q^ A for zth process are 
determined by the instantaneous charging states of two 
relevant granules and by the local electric field F n and 
temperature T, accordingly to the expressions: 

?£i = (1 - <£) (1 - <£+a) f (eF n • A + E c ) 

Ini = ^(l-^ +A )^(-^nF n -A) 



9n,A = 2 CTnCrn + A (^n^n+A - 1) X 

x ip (eer n+A F n ■ A - E c ) . 



(2) 



Thus the charging energy is positive, E c , for the pair 
creation, zero for the transport, and negative, — E c , 
for the recombination processes. The function tp(E) = 
ujNpE/[cxp(f3E) — 1] expresses the total probability, at 
given inverse temperature (i = l/(fceT), for electron 
transition between granules with Fermi density of states 
N-p and Fermi levels differing by E. The hopping fre- 
quency lo — uj a exp(— 2\s) involves the attempt frequency, 
uj a ~ Ep/H, the inverse tunneling length \ (typically 
~ 10 urn -1 ), and the inter-granule spacing s = a — d. 
Local electric field F n on nth site consists of the external 
applied field A (site independent) and the Coulomb field 
C n due to all other charges in the system: 



C n 



£ eff 



n' 



n 



(3) 



A suitable approximation is achieved with passing from 
discrete-valued functions cr n of discrete argument n = 
a(ni,n 2 ) to their continuous- valued mean-field (MF) 
equivalents a r = (a n ) r (mean charge density) and p r = 
(<J n ) r (mean charge carrier density). These densities are 
obtained by averaging over a wide enough area (that is, 
great compared to the lattice period but small compared 
to the size of entire system or its parts) around any point 
r in the plane (for simplicity, we drop the position index 
at averages () r in what follows). This also implies pass- 
ing to a smooth local field: 



A+- 



e eff a' 



r 



- r|3 



dr'. 



(4) 



and to the averaged transition rates q^\ — (q^ A ) and 



p£a = ^"Ca^- These rates fully define the temporal 
derivatives of mean densities: 



) __(2) , (2) 

+A,-A ^r.A ^*V+A,-A 



D (3) 
Pr.A 



(5) 



V \a W +o (1) -a (2) +a (2) 

/ / [y r ,A + y r +A,-A y r ,A + 1t+a,-a 



III. MEAN-FIELD DENSITIES IN 
EQUILIBRIUM 

We perform the above defined averages in the simplest 
assumption of no correlations between different sites: 
(/nflw) = (fn) (aw), n' ^ n, and using the evident rules: 
(<J n k+1 ) = cr r , (<J n k ) = p r - The resulting averaged rates 
are: 



o (1) 
o (2) 

1r,A 



(2) 
Pr,A 



o (3) 

1r,A 



„(3) 
Pr.A 



a° r a° r+A p(eF r -A + E c ), 
a° r+A [a+<p(-eF T - A) 

+ <7-y>(eF P -A)] , 

a r ° +A [a+tp (-eFr- A) 

- a~ip (eF r • A)] , 
[cr+<j; +AV (-eF r -A-E c ) 

+ a-a+ +A ip (eF r • A - E c )] , 
[a+a- +A ip(-eF r -A-E c ) 

-v7Vr + A<P(eFr-A-E c )], (7) 



where the mean occupation numbers for each charging 
state af = (p r ± cr r )/2 and of = 1 — p r satisfy the nor- 
malization condition: Y] ■ a % r = 1. 

In a similar way to Eq. 5, we express the vector of 
average current density j n at nth site: 



+A,-A 



+ D (2) -D (2) +» (3) 



(8) 



and then its MF extension j r is obtained by simple re- 
placing n by r in the arguments of qrM and pW. Ex- 
panding these continuous functions in powers of |A| = a, 
we conclude that Eq. 5 gets reduced to usual continuity 
equation: 



Or = V 2 • Jr, 

e 



(9) 



with the two dimensional (2D) nabla: V 2 = (d x ,d y ). 
We begin the analysis of Eqs. 5-9 from the simplest 
situation of thermal equilibrium in absence of electric 
field, F r = 0, then Eq. 5 turns into evident identity: 
(T r = 0, that means zero charge density, and Eq. 8 yields 
in zero current density: j r = 0, while Eq. 6 provides a 
finite and constant value of charge carrier density: 



(10) 



2 + exp (pE c /2) ' 



9r,A 



(6) 



The set of Eqs. 2-6 provides a continuous description of 
the considered system, once a proper averaging procedure 
is established. 



At low temperatures, f3E c 3> 1, this value is exponen- 
tially small: p e rts 2cxp(—(3E c /2), and for high tempera- 
tures, j3E c <C 1, it behaves as p e w p^ — (3E c /9, tending 
to the limit poo = 2/3, corresponding to equipartition 
between all three fractions a % (Fig. 4, though this limit 
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FIG. 4: Equilibrium density p e of charge carriers in func- 
tion of temperature (solid line). The curve 1 (dashed 
line) corresponds to the low temperature asymptotics p e ?a 
2exp (— E c /2k-g.T), and the curve 2 (dash-dotted line) to the 
high temperature asymptotic p e « pao — E c /9kBT, converging 
to the limit poo = 2/3 (dotted line). 
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FIG. 5: The charge density a in function of the carrier density 
p for different temperatures (corresponding to different ther- 
mal equilibrium values p e ). Note closeness of all the curves to 
that for low-temperature limit p e — > 0, given exactly by Eq. 
12. 



being beyond actual validity of the model, as indicated 
in Sec. II). 

In presence of electric fields F r ^ 0, the local equilib- 
rium should be perturbed and the system should generate 
current and generally accumulate charge. Then, from Eq. 
6, the charge density a r is related to the carrier density 
p r as: 



2 {Pr - Pe) (Pr + Pe - 2pePr) 
<J r = 



(1-PeY 



(11) 



describing the increase of charge density with going away 
from equilibrium. As seen from Fig. 5, for not too high 
temperatures T < E c /kB where the neglect of multiple 
charged states is justified, this dependence is reasonably 
close to the simplest low-temperature form: 



V? 



Pf: 



(12) 



that will be practically used in what follows. 

Now we are in position to pass to the out-of- 
equilibrium situations, beginning from a simpler case of 
dc current flowing through the FA. 



IV. STEADY STATE CONDUCTION IN FA 

In presence of (generally non-uniform) fields F r and 
densities <7 r , p r , we expand Eq. 8 up to 1st order terms 
in |A| = a and obtain the local current density as a sum 
of two contributions, the field-driven and diffusive: 



t ld +if =9(Pr)F r 



eD(p T )V 2 a T , (13) 



where the effective conductivity g and diffusion coefficient 
D are functions of the local charge carrier density, p= p r : 



g{p) = - 2(l-p)V(£c)+/>(l-pV(0) 
1 



^(p 2 -aV(-£ C ) 



p(l - p e )V(0)(l - p)p 2 M-Ec)/2 



p(l - 2 Pe ) + pi 



(14) 



In view of Eqs. 11, 12, we can consider g and D as 
even functions of local charge density a, and just this 
dependence will be mostly used below. Also g and D 
depend on temperature, through the functions ip and if' . 
The system of Eqs. 11 -14, together with Eq. 4, is closed 
and self-consistent, defining the distributions of a r and p r 
at given j r . It is readily seen to admit the trivial solution, 
cr(x) = 0, and now we shall argue that in fact this is the 
only practical solution for FA. 

First of all, we notice physical restrictions on the 
charge accumulation in FA. By the problem symmetry, 
the charge density should only depend on the coordi- 
nate along the current, a = cr(x), this function be- 
ing odd (in the geometry of Fig. 3) and supposedly 
monotonous. Then its maximum value cr max = cx(L/2) 
will define the characteristic scale for the Coulomb field: 
C <~ cr max e/ (e c fla 2 ) which should not be higher than 
typical applied fields A ~ 10 2 V/cm (as seen from rela- 
tively moderate non-ohmic vs ohmic response in the ex- 
periment). Thus the maximum charge density should not 
surpass the level of Ae e ga 2 /e <~ 10~ 3 , that is much lower 
than the equilibrium density of charge carriers p e (except 
for, maybe, too low temperatures, T < 0.07E c /kQ <~ 10 
K). Therefore, one can neglect the small difference, Eq. 
12, setting constant values: p m p e and then g w g e = 
g( Pe ),D^D e = D( Pe ). 



Under such condition, we can eliminate the (not well 
known) constant A from Eq. 13, bringing this equation 
to the integro-differential form: 



d 2 a(x) 
dx 2 



<Je 



D e e eS a 2 



P 



1/2 o{x')dx' 



-1/2 



r /\2 ■ 



(15) 



where the P-symbol at integration in x' means the "dis- 
crete principal value", that is omission of the interval 
[x — a. x + a) to avoid the apparent divergence, in agree- 
ment with the minimum distance between granules in the 
lattice. Thus the regularized integral converges rapidly, 
then it is reasonable to fix the argument of a-density at 
x' = x, arriving at a simple differential equation: 



d 2 a(x) _ a(x) 
dx 2 r 2 , 



(16) 



Here the parameter 



4 



a 3 2eP E ° + 5e^/ 2 + 2 
d e 3/3s c /2 + 2pE c <eP E <> - e^/ 2 



defines the temperature dependent length scale rp, and 
the x-odd solution of Eq. 16 is just o~(x) = o\ smh(x/rp). 
However, for all the considered temperatures, (3E C > 1 
(see the note in Sec. II), this scale is rp < a, that is 
by many orders of magnitude smaller than the FA size I. 
Then the estimate for the constant G\ in the above solu- 
tion, (Ti ~ (J max e~ l / rf> with the exponent as great as for 
instance l/rp ~ 10 4 , makes this solution practically van- 
ishing within whole FA, except maybe for a very narrow 
vicinity ~ rp of its interface with CA (where, strictly 
speaking, Eq. 16 no more holds). This evident conse- 
quence of long-range character of Coulomb fields in FA 
will be contrasted below with the situation in CA, where 
charge accumulation turns possible due to screening ef- 
fects by the metallic contacts and to the related short- 
range fields. 

Thus we conclude that there is practically no charge 
accumulation and hence no diffusive contribution to the 
current in FA. Thus the steady state of FA in out-of- 
equilibrium conditions should be characterized by the 
ohmic conductivity g e . In fact, an estimation (based on 
an experimental system [17]) suggests that the FA con- 
tribution to the overall resistance turns to be about two 
orders of magnitude smaller than the CA one (see be- 
low), and thus the transport is expected to be mainly 
controlled by CA. 



V. STEADY STATE CONDUCTION IN CA 
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FIG. 6: Kinetic processes between nth granule and the metal- 
lic electrode in CA. 



granule and, using the same techniques that before, their 
mean values are: 



ll 4) = 

<$» = 

& = 

J 7 ) - 



(Pr ~ <T,) 1p(U ~ E' c ) , 

(l-p T )i>(U+E' c ), 
(1- Pr )i>(-U + E' c ). 



(17) 



Here the function ip(E) formally differs from tp(E) only 
by changing the pre-factor: uj — > lu 1 = u) a e~ 2xb <C u), but 
the arguments of these functions in Eq. 17 include other 
characteristic energies. Thus, the energy U — eb'S is due 
to the electric field S = F c (z = b') at the contact surface 
above the granule. As seen from Fig. 7, this field is al- 
ways normal to the surface and its value is defined by the 
local charge density a (see below). At least, the charging 
energy E' c for a granule under the contact can be some- 
what lower (e.g., by ~ 1/2) than E c . Then the kinetic 
equations in interface region present a generalization of 
Eqs. 5-6, as follows: 



y r +A,-A 



7( 4 ) 



7( 5 ) 



J 2 ) 
Pr.A 



,(6) 



„(2) 

Pr+A,-A 



,(7) 



D (3) 
Pr,A 



(18) 



V \a (1) -o (2) +o (2) -o (3) 

Z / L^r,A + ?r+A,-A 1r,A + «r+A,-A ?r,A 



A 



,(4)-„(5) +g C6) +? (7) 



(19) 



The additional terms, by the normal processes 4) to 7), 
are responsible for appearance of a normal component of 
current density: 



The kinetics in CA includes, besides the processes 1) 
to 3) of Sec. Ill and IV, also four additional microscopic 
processes between nth granule and the electrode (Fig. 6) 
which are just responsible for variations of total charge 
Q by ±1. The respective rates i = 4, . . .7, are also 
dependent on the charging state (cr r ,p r ) of the relevant 



7(4) 



(5) 



,(6) 



7(7) 



(20) 



besides the planar component, still given by Eq. 8. But 
an even more important difference from the FA case is 
the fact that the Coulomb field here is formed by a dou- 
ble layer of charges, those by granules themselves and by 



FIG. 7: Formation of local electrical fields by a dipole of a 
charged granule and its (oppositely charged) image: at the 
surface of the metallic electrode (point a) and on other gran- 
ule (point 6). 



their images in the metallic electrode (Fig. 7) . Summing 
the contributions from all the charged granules and their 
images (except for the image of nth granule itself, al- 
ready included in the energy E' c ), we find that the above 
mentioned field at the contact surface above the point 
r of the granular layer, S^, can be expressed as a local 
function of the charge density <r r : 



S r = C r {z = b') 



47TC 



(21) 



replacing the integral relations, Eqs. 3-4, in FA. Also, 
note that the relevant dielectric constant for this field 
formed outside the granular layer is rather the host value 
e than the renormalized e e f / within the layer (as by Eq. 
3). Then, the planar component of the field by charged 
granules F£ ( = C r (z — 0) is determined by the above de- 
fined normal field S r through the relation = &'V2-5 r - 
The density of planar current is j£ z = gF? 1 — e£>V2fr, 
accordingly to Eq. 13, that is both field-driven and dif- 
fusive contributions into j^' are present here and both 
they are proportional to the gradient of a r . In the low 
temperature limit, this proportionality is given by: 



.pi 



ine 3 ujN F b' , , euN F k B T 
^—9 Or) + 



£ e //cr a 
Note that the presence of a non-linear function: 



(22) 



50) = "JpI + • 



2 P l 



defines a non-ohmic conduction in CA. In fact, this func- 
tion should be defined by Eq. 21 only for charge density 
below its maximum possible value |cr„ mx | = \J 1 — p 2 , 
turning zero for \a\ > \(r max \ (note that the latter re- 
striction just corresponds to our initial limitation to 
the single charged states, see Sec. II). In the same 
limit of low temperatures, the normal current density 
is obtained from Eqs. 16, 17 as j z (r) = G z "E r where 
G z sa (jj'N-pE' c E e ff /An. Finally, the kinetic equation in 
this case is obtained, in analogy with Eq. 8, as: 



a 2 b. 



-V a j 



(23) 



FIG. 8: Relations between longitudinal (j x ) and normal (j z 
currents in CA, adding to the total current I. 



This equation permits to describe the steady state con- 
duction as well as various time dependent processes. The 
first important conclusion is that steady state conduction 
in the interface turns only possible at non-zero charge 
density gradient, that is, necessarily involving charge ac- 
cumulation, in contrast to the above considered situation 
in bulk. 

Let us restrict here the analysis to the steady state 
conduction regime which is simpler, though the obtained 
results can be also used for the analysis of a more involved 
case when an explicit temporal dependence of charge den- 
sity is included in Eq. 23 (this will be a topic of future 
study). 

We choose the contacts geometry in the form of a rect- 
angular stripe of planar dimensions L x L', along and 
across the current respectively. In neglect of relatively 
small effects of current non-uniformity along the lateral 
boundaries, the only relevant coordinate for the problem 
is longitudinal, x (Fig. 8), so we consider the relevant 
function a x with its derivatives, spatial o~ x and temporal 
d x . In the steady state regime, & = in Eq. 23, and the 
total current I = const, defined by the action of external 
source. Then, using the above approximation for g(a), a 
non-linear 2nd order equation for charge density is found: 



dx 



{[9 (O + rR} 



k 2 c 



0. 



(24) 



Here the parameters are: k 2 = (us'E' c )/ (abLok-oTi) and 
t = T/Ti, where T is the actual temperature and T\ = 
8ire 2 b' /a 2 k-Q£ef f- To define completely its solution, the 
following boundary conditions are imposed: 



'x=0 



k 2 b'a x=0 
g(<J x =o) +t' 



I 



LeujbN F k B Ti g (a x = L ) + r 



(25) 



(26) 



Here Eq. 25 corresponds to the fact that the longitu- 
dinal current j x at the initial point of contact/granular 
sample interface (the leftmost in Fig. 8) is fully supplied 
by the normal current j z entering from the contact to the 
granular sample, and Eq. 26 corresponds to the current 
continuity at passage from CA (of length L along the x 
axis) to FA. 
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Let us discuss the solution of Eq. 24 qualitatively. 
Generally, to fulfill the conditions, Eqs. 25, 26, one 
needs a quite subtle balance to be maintained between 
the charge density and its derivatives at both ends of con- 
tact interface. But the situation is radically simplified 
when the length L is much greater than the characteris- 
tic decay length for charge and current density: kL >• 1. 
In this case, the relevant coordinate is £ = L — x, so that 
the boundary condition 25 corresponds to £ = L — > oo, 
when both its left and right hand side turn zeros: 



0. 



0. 



(27) 



The numeric solution shows that, for any initial (with 
respect to £, that is related to x = L, Eq. 26) value of 
charge density cre = o = <jq, there is a unique initial value 
of its derivative cr'^_ = D(a ) which just assures the lim- 
its, Eq. 27, while for cr^_ > D{uq) the asymptotic value 
diverges as er^oo — > oo, and for cr^ =0 < D(a ) it diverges 
as cr^oo — > — oo. Then, using the boundary condition, 
Eq. 26, and taking into account the relation V = Voao 
following from Eq. 23 with Vb = 47refe'/(£ e //a 2 ), we con- 
clude that the function D(a ) generates the I-V charac- 
teristics: 



V 



(28) 



where I\ = ecoNpk^Ti. 

A more detailed analysis of Eq. 24 is presented in 
Appendix. In particular, for the weak cu rrent regime 
(Regime I) when ctq <C <J\ = y / 32p e (p e + t) < 1. so that 
g[a) ?a p e + a 2 / (2p e ) along whole the contact interfaces, 
Eq. 23 admits an approximate analytic solution: 



<7£ = cr e 



-2A£ 



(29) 



with the exponential decay index A = k/y/ p e + r. 

This results in the explicit I-V characteristics for 
Regime I: 



I = G V 



(30) 



for V < Vi = <J\Vq, Eq. 30 describes the initial ohmic 
CA conductance (temperature r dependent): 



hkb' 

Go = -^—VPe {r)+r, 



(31) 



which turns non-ohmic for V ~ V\ . But at so high volt- 
ages another conduction regime already applies (called 
Regime II), where u\ <C <Jq <C 1 and one has g(a) w a 
(see Eq. 21). Following the same reasoning as for the 
Regime I, we obtain a non-linear I- V characteristics for 
Regime II: 



(32) 



this law is weaker temperature dependent than Eq. 30, 
which is related to the fact that the conductance in 
Regime II is mainly due to dynamical accumulation of 
charge and not to thermic excitation of charge carriers. 
Interestingly a I cx V 3 ^ 2 law was recently found in exper- 
imental measurements [17] . Further, such non-linearity 
can be yet more pronounced if multiple charging states 
are engaged, as may be the case in real granular lay- 
ers with a certain statistical distribution of granule sizes 
present. 

At least, for even stronger currents, when already cr <~ 
1, the solutions of Eq. 24 can be obtained numerically, 
following the above discussed procedure of adjustment of 
the derivative D(ao) to a given ctq- Such solutions have 
an asymptotic behavior of the type: I cx V 5 / 4 . 

A simple and important exact relation for the total 
accumulated charge Q in CA is obtained from the direct 
integration of Eq. 24: 

Q = tl, 

where the parameter t = l/tp(—E' c ) should have a role 
of characteristic relaxation time in non-stationary pro- 
cesses. Assuming its value t~ls (comparable with the 
experimental observations [11]), together with the above 
used values of ui and T\ , we conclude that the character- 
istic length scale A -1 for solutions of Eq. 24 can reach 
up to <~ 10 3 a ~ 1/im, which is a reasonable scale for a 
charge distribution beneath the contacts. 



VI. GLOBAL CONDUCTION IN THE SYSTEM 

The conduction in the overall system results from 
matching of the above considered processes in CA and 
FA. Thus, in order to evaluate the global resistance of 
this circuit in series it is necessary to add the contri- 
butions of both areas to it. Recent measurements [17] 
have shown notably non-linear I-V curves (already at low 
enough voltages), so, accordingly to the above discussion, 
this indicates that the resistance should be dominated by 
CA. To have a clear view on it, we can use the typical 
parameters for the granular film: a <~ 5 nm, d <~ 4 nm, 
X ~ 10 nm , 6^8 nm, b' <~ 2 nm, E c ~ 10 meV, 
Np <~ 1 eV _1 and take w as a (less known) fitting pa- 
rameter. For the considered rectangular CIP geometry 
we also use the experimental values [17] of width L' = 3 
mm and of distance between the contacts I = 100 pm. 

Choosing T = 50 K, the ohmic conductance of the 
FA, Gpa, can be calculated through the formula Gpa — 
g(p e )L'b/l ~wl.5x 10~ 18 S. In the CA, we can estimate 
the conductance (in Regime I) following the above for- 
mula Gca ~ w 8.0 x 10~ 22 S. Thus it is clear that, for any 
choice of uj, the conductance of the CA is about 4 orders 
of magnitude smaller than that of the FA and for that 
reason it should dominate the global resistance of the sys- 
tem. Then, using the formulae, Eqs. 29-31, we obtain a 
good agreement with the experimental data by Ref. [17] 
as shown in Fig. 9. It should be noted however that the 
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FIG. 9: I-V characteristics for a granular sample at different 
temperatures, compared with the theoretical curves for for 
regimes I and II. Inset: temperature dependence of ohmic 
conductance Go, measured data (circles) vs calculated by Eq. 
31. 



accumulation, leading to a non-ohmic / oc V 3//2 conduc- 
tion law. The calculated I-V curves and temperature 
dependencies are found in a good agreement with avail- 
able experimental data. The proposed model can be fur- 
ther developed for description of multilayer strucuture ef- 
fects and also of non-stationary conduction processes, like 
anomalous slow current relaxation [13]. Finally, the elas- 
tic effects of Coulomb forces by charged granules can be 
included in order to explain the remarkable phenomenon 
of resistive-capacitive switching [18], in granular layered 
conductors. 
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effective value of the parameter Vo giving the best fit to 
the experimental data should be notably higher then that 
given by our formula (before Eq. 28) for single layer sys- 
tem. Thus, with the above choice of other parameters, we 
have the single-layer value Vq ~ 0.5 V whereas the best fit 
for 10-layer experimental sample needs instead V exp rts 3 
V. This difference can be effectively accounted for by a 
simple multiplicative factor a w 6 (the "multilayer fac- 
tor") so that V exp = aVo assures both the agreement for 
Regimes I, II of I-V curves and the boundary V <~ V exp 
between them, clearly seen in Fig 9. 



VII. CONCLUSION 

In conclusion, the mean-field model is developed for 
tunnel conduction in a granular layer, including three 
principal processes of creation and annihilation of pairs 
of opposite charges on neighbor granules and of charge 
transfer from a charged granule to a neighbor neutral 
granule. Effective kinetic equations for averaged charge 
densities are derived for the characteristic areas of the 
granular sample: the contact areas beneath metallic cur- 
rent leads and free area between these leads. From these 
kinetic equations, it is shown that the tunnel conduc- 
tion in the free area does not produce any notable charge 
accumulation, and the conduction regime here is purely 
ohmic. Contrariwise, such conduction in the contact area 
turns impossible without charge accumulation, leading to 
generally non-ohmic conduction regime, since the contact 
area dominates in the overall resistance. Approximate 
analytic treatment is developed for calculation of charge 
density and tunnel current in two characteristic regimes: 
I) for weak charge accumulation (compared to the ther- 
mal density of charge carriers) and II) for strong charge 



IX. APPENDIX 

Let us consider the equation: 

Ti [g{a) + T] % - fcV = ( A1 ) 

with certain boundary conditions tr(0) = a , ct'(0) = a' , 
resulting from Eqs. 24, 25. For a rather general function 
g(a) we can define the function 

f(a)= f g(a')da', (A2) 
Jo 

then Eq. Al presents itself as: 

^ = fc2 ^> ( A3 ) 

where = f (cr^) + rcr^. Considered irrespectively of £: 

f(a)+Ta = F, (A4) 

this equation also defines a as a certain function of F: 
a = a(F). Hence it is possible to construct the following 
function: 

4>{F) = 2 f (j(F')dF'. (A5) 
Jo 

Now, multiplying Eq. A3 by 2c?F/d£, we arrive at the 
equation: 

dn\dn) - k de' (A6) 
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FIG. 10: Charge density and current distribution in the CA 
region (Regime I). 




FIG. 11: Charge density distribution in Regime II. A fast 
decay is changed to a slower exponential law, after density 
dropping below the characteristic value p e . 



with 0(£) = 4>(F^). Integrating Eq. A6 in £, we obtain a 
1st order separable equation for F^: 



(A7) 



We expect the function F to decrease at going from 
£ = into depth of interface region, hence choose the 
negative sign on r.h.s. of Eq. A7 and obtain its explicit 
solution as: 



rFc, 
J F e 



dF' 



VW 7 



= kii 



(A8) 



with F = / ((To) + t<t . Finally, the sought solution 
for <j£ = a (F^) results from substitution of the function 
F ? , given implicitly by Eq. A8, into <r(F) defined by 
Eq. A4. Consider some particular realizations of the 
above scheme. For the approximate solution of g(u) given 
above, we have the explicit integral, Eq. A2, in the form: 



F(a) = f(a)+ra 



T + 



(A9) 



In the case a <C p e <C 1 (Regime I), Eq. A9 is approxi- 
mated as: 



F w (p e + r) a 



6p e 



(A10) 



hence <r(F) corresponds to a real root of the cubic equa- 
tion, Eq. A10, and in the same approximation of Regime 



I it is given by: 

a(F) 



F 



l-^r), (Ail) 



with (7i = 4\J p e (p e + t) 3 - Using this form in Eq. A5, 
we obtain: 



<p{F) 



Pe +T 



4F J 

1 - "3- ) - (A12) 



and then substituting into Eq. A8: 



ln-t 



l + ^/l-(2F/(T 1 ) 2 



F 



1 - {2F /a 1 y 



(A13) 



Inverting this relation, we define an explicit solution for 



F(0 « F c- A « 



1 + 3(1- e" a *) 



(A14) 



Finally, substituting Eq. A14 into Eq. All, we arrive at 
the result of Eq. 29 corresponding to Fig. 10. 
For the regime II we have in a similar way: 

F(a) « <t(t + */2), 
a{F) M y/2F + t 2 - r, 



(2F + r 2 ) 3/2 -r(3F + r 2 ) 



,1/4 



3t 



2V4 (F 1/4 -A^ 



(A15) 
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with Ai = fc/(2 3 / 4 v / 3), obtaining the charge density dis- 
tribution (Fig. 11): 

«(0 ~ (V^+T - XiS,) 2 - r. (A16) 
This function seems to turn zero already at £ = 



{\f a o + T ~~ V^)/^ii °ut hi f ac t the fast parabolic decay 
by Eq. A16 only extends to £ <~ £*, such that cr^» ~ p e , 
and for £ > £* the decay turns exponential, like Eq. 29. 
The /-V characteristics, Eq. 32, follows directly from 
Eq. A16. 
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